{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pyamg.gallery import diffusion_stencil_2d, stencil_grid\n",
    "from pyamg.classical.split import RS, PMIS, MIS\n",
    "from scipy.io import mmwrite\n",
    "from scipy.sparse import csr_matrix\n",
    "import PetscBinaryIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PetscPrint(filename, Mat):\n",
    "    PrintMat = csr_matrix(Mat)\n",
    "    PetscBinaryIO.PetscBinaryIO().writeMatSciPy(open(filename,'w'), PrintMat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ClassicStrength(A, theta = 0.0):\n",
    "    A = A.tocsr()\n",
    "    \n",
    "    Sp = np.empty_like(A.indptr)\n",
    "    Sj = np.empty_like(A.indices)\n",
    "    Sx = np.empty_like(A.data)\n",
    "    \n",
    "    nnz = 0\n",
    "    n = A.shape[0]\n",
    "    \n",
    "    Sp[0] = 0\n",
    "    for i in range(0, n):\n",
    "        diag = 0\n",
    "        for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "            if (A.indices[j] == i):\n",
    "                diag = A.data[j]\n",
    "                break\n",
    "        \n",
    "        if diag < 0:\n",
    "            row_scale = -np.inf\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if (val > row_scale):\n",
    "                    row_scale = val\n",
    "        else:\n",
    "            row_scale = np.inf\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if (val < row_scale):\n",
    "                    row_scale = val\n",
    "        \n",
    "        threshold = row_scale * theta\n",
    "        \n",
    "        Sj[nnz] = i\n",
    "        Sx[nnz] = diag\n",
    "        nnz += 1\n",
    "        \n",
    "        if diag < 0:\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if val > threshold:\n",
    "                    col = A.indices[j]\n",
    "                    Sj[nnz] = col\n",
    "                    Sx[nnz] = val\n",
    "                    nnz += 1\n",
    "        else:\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if val < threshold:\n",
    "                    col = A.indices[j]\n",
    "                    Sj[nnz] = col\n",
    "                    Sx[nnz] = val\n",
    "                    nnz += 1\n",
    "                    \n",
    "        Sp[i + 1] = nnz\n",
    "\n",
    "    return csr_matrix((Sx, Sj, Sp), shape=A.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SymmetricStrength(A, theta = 0.0):\n",
    "    A = A.tocsr()\n",
    "    \n",
    "    Sp = np.empty_like(A.indptr)\n",
    "    Sj = np.empty_like(A.indices)\n",
    "    Sx = np.empty_like(A.data)\n",
    "    \n",
    "    n = A.shape[0]\n",
    "    \n",
    "    diags = np.zeros(n)\n",
    "    row_thresholds = np.zeros(n)\n",
    "\n",
    "    nnz = 0\n",
    "    Sp[0] = 0\n",
    "    for i in range(0, n):\n",
    "        diag = 0\n",
    "        for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "            if (A.indices[j] == i):\n",
    "                diags[i] = A.data[j]\n",
    "                break\n",
    "        \n",
    "        if diag < 0:\n",
    "            row_scale = -np.inf\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if (val > row_scale):\n",
    "                    row_scale = val\n",
    "        else:\n",
    "            row_scale = np.inf\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                if (val < row_scale):\n",
    "                    row_scale = val\n",
    "        \n",
    "        row_thresholds[i] = row_scale * theta\n",
    "    \n",
    "    for i in range(0, n):\n",
    "        Sj[nnz] = i\n",
    "        Sx[nnz] = diags[i]\n",
    "        nnz += 1\n",
    "        \n",
    "        if diag < 0:\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                col = A.indices[j]\n",
    "                if val > row_thresholds[i] or val > row_thresholds[col]:\n",
    "                    Sj[nnz] = col\n",
    "                    Sx[nnz] = val\n",
    "                    nnz += 1\n",
    "        else:\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                val = A.data[j]\n",
    "                col = A.indices[j]\n",
    "                if val < row_thresholds[i] or val < row_thresholds[col]:\n",
    "                    Sj[nnz] = col\n",
    "                    Sx[nnz] = val\n",
    "                    nnz += 1\n",
    "                    \n",
    "        Sp[i + 1] = nnz\n",
    "\n",
    "    return csr_matrix((Sx, Sj, Sp), shape=A.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ModClassicalInterp(A, S, splitting):\n",
    "    A = A.tocsr()\n",
    "    S = S.copy()\n",
    "    S = S.tocsr()\n",
    "\n",
    "    Pp = np.empty_like(A.indptr)\n",
    "    \n",
    "    col_to_new = np.empty(A.shape[0])\n",
    "    col_to_new.fill(-1)\n",
    "\n",
    "    nnz = 0\n",
    "    Pp[0] = nnz\n",
    "    for i in range(A.shape[0]):\n",
    "        if splitting[i] == 1:\n",
    "            nnz += 1\n",
    "        else:\n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                if splitting[col] == 1:\n",
    "                    nnz += 1\n",
    "        Pp[i+1] = nnz\n",
    "                \n",
    "    Pj = np.empty(nnz, dtype=Pp.dtype)\n",
    "    Px = np.empty(nnz, dtype=A.dtype)\n",
    "    \n",
    "    nnz = 0\n",
    "    ctr = 0\n",
    "    \n",
    "    A.sort_indices()\n",
    "    S.sort_indices()\n",
    "    \n",
    "    row_strong = np.zeros(A.shape[0], dtype='float')\n",
    "    pos = np.empty(A.shape[0], dtype='int')\n",
    "    pos.fill(-1)\n",
    "            \n",
    "    n_cols = 0\n",
    "    for i in range(A.shape[0]):\n",
    "        if splitting[i] == 1:\n",
    "            Pj[nnz] = i\n",
    "            Px[nnz] = 1.0\n",
    "            col_to_new[i] = n_cols\n",
    "            n_cols += 1\n",
    "            nnz += 1\n",
    "        else:            \n",
    "            # Add coarse points in row\n",
    "            ctr = S.indptr[i]\n",
    "            endS = S.indptr[i+1]\n",
    "            sign = 1\n",
    "            weak_sum = 0\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                col = A.indices[j]\n",
    "                val = A.data[j]\n",
    "                if ctr < endS and S.indices[ctr] == col:\n",
    "                    if splitting[col] == 1:\n",
    "                        pos[col] = nnz\n",
    "                        Pj[nnz] = col\n",
    "                        Px[nnz] = val\n",
    "                        nnz += 1\n",
    "                    elif col == i:\n",
    "                        weak_sum += val\n",
    "                        if val < 0:\n",
    "                            sign = -1\n",
    "                    else:\n",
    "                        row_strong[col] = val\n",
    "                    ctr += 1\n",
    "                else:\n",
    "                    weak_sum += val\n",
    "    \n",
    "            # Add distance two to coarse sum \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                if col == i or splitting[col] == 1:\n",
    "                    continue\n",
    "                coarse_sum = 0\n",
    "                for k in range(A.indptr[col], A.indptr[col+1]):\n",
    "                    colk = A.indices[k]\n",
    "                    if pos[colk] >= 0:\n",
    "                        val = A.data[k]\n",
    "                        if val * sign < 0:\n",
    "                            coarse_sum += val\n",
    "                if coarse_sum == 0:\n",
    "                    weak_sum += S.data[j]\n",
    "                    row_strong[col] = 0\n",
    "                else:\n",
    "                    row_strong[col] /= coarse_sum\n",
    "            \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                row_val = row_strong[col]\n",
    "                if col == i or splitting[col] == 1:\n",
    "                    continue\n",
    "                \n",
    "                if row_val != 0:\n",
    "                    for k in range(A.indptr[col], A.indptr[col+1]):\n",
    "                        colk = A.indices[k]\n",
    "                        idx = pos[colk]\n",
    "                        if idx >= 0:\n",
    "                            val = A.data[k]\n",
    "                            if val * sign < 0:\n",
    "                                Px[idx] += row_val * val \n",
    "                                \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                idx = pos[col]\n",
    "                if idx >= 0:\n",
    "                    Px[idx] /= -weak_sum\n",
    "                    pos[col] = -1\n",
    "                else:\n",
    "                    row_strong[col] = 0\n",
    "                    \n",
    "    for i in range(nnz):\n",
    "        Pj[i] = col_to_new[Pj[i]]\n",
    "                    \n",
    "    return csr_matrix((Px, Pj, Pp), shape = (A.shape[0], n_cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ExtendInterpolation(A, S, splitting):\n",
    "    A = A.tocsr()\n",
    "    S = S.copy()\n",
    "    S = S.tocsr()\n",
    "    A.sort_indices()\n",
    "    S.sort_indices()\n",
    "    \n",
    "    Pp = np.empty_like(A.indptr)\n",
    "    pos = np.empty(A.shape[0], dtype='int')\n",
    "    pos.fill(-1)\n",
    "    \n",
    "    nnz = 0\n",
    "    Pp[0] = nnz\n",
    "    for i in range(A.shape[0]):\n",
    "        if splitting[i] == 1:\n",
    "            nnz += 1\n",
    "        else:\n",
    "            row_pos = list()\n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                if col == i: \n",
    "                    continue\n",
    "                \n",
    "                # Add distance one and two connections\n",
    "                if splitting[col] == 1:\n",
    "                    if pos[col] == -1:\n",
    "                        pos[col] = 1\n",
    "                        row_pos.append(col)\n",
    "                        nnz += 1\n",
    "                else:\n",
    "                    for k in range(S.indptr[col], S.indptr[col+1]):\n",
    "                        colk = S.indices[k]\n",
    "                        if splitting[colk] == 1 and pos[colk] == -1:\n",
    "                            pos[colk] = 1\n",
    "                            row_pos.append(colk)\n",
    "                            nnz += 1\n",
    "            for col in row_pos:\n",
    "                pos[col] = -1\n",
    "        Pp[i+1] = nnz\n",
    "                \n",
    "    Pj = np.empty(nnz, dtype=Pp.dtype)\n",
    "    Px = np.empty(nnz, dtype=A.dtype)\n",
    "    \n",
    "    nnz = 0\n",
    "    ctr = 0\n",
    "    \n",
    "    row_strong = np.zeros(A.shape[0], dtype='float')\n",
    "            \n",
    "    for i in range(A.shape[0]):\n",
    "        if splitting[i] == 1:\n",
    "            Pj[nnz] = i\n",
    "            Px[nnz] = 1.0\n",
    "            nnz += 1\n",
    "        else:            \n",
    "            # Add coarse points in row\n",
    "            weak_sum = 0\n",
    "            \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                val = S.data[j]   \n",
    "                \n",
    "                if col == i:\n",
    "                    continue\n",
    "                    \n",
    "                # Add Distance-One Connections\n",
    "                if splitting[col] == 1: \n",
    "                    if pos[col] == -1:\n",
    "                        pos[col] = nnz\n",
    "                        Pj[nnz] = col\n",
    "                        Px[nnz] = val\n",
    "                        nnz += 1\n",
    "                    else:\n",
    "                        Px[pos[col]] = val    \n",
    "                else:\n",
    "                    row_strong[col] = val\n",
    "                        \n",
    "                    # Add distance two connections\n",
    "                    for k in range(S.indptr[col], S.indptr[col+1]):\n",
    "                        colk = S.indices[k]\n",
    "                        if splitting[colk] == 1 and pos[colk] == -1:\n",
    "                            pos[colk] = nnz\n",
    "                            Pj[nnz] = colk\n",
    "                            Px[nnz] = 0\n",
    "                            nnz += 1\n",
    "            \n",
    "            # Add weak connections not in C^_{i} to weak_sum\n",
    "            weak_sum = 0\n",
    "            ctr = S.indptr[i]\n",
    "            endS = S.indptr[i+1]\n",
    "            for j in range(A.indptr[i], A.indptr[i+1]):\n",
    "                col = A.indices[j]\n",
    "                val = A.data[j]  \n",
    "                if ctr < endS and S.indices[ctr] == col:\n",
    "                    ctr += 1\n",
    "                    if col == i:\n",
    "                        weak_sum += val\n",
    "                else:\n",
    "                    # Weak sum\n",
    "                    if pos[col] == -1:\n",
    "                        weak_sum += val\n",
    "            \n",
    "            # Add distance two to coarse sum \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                \n",
    "                # Find diag val\n",
    "                for k in range(S.indptr[col], S.indptr[col+1]):\n",
    "                    if S.indices[k] == col:\n",
    "                        val = S.data[k]\n",
    "                        break\n",
    "                sign = 1\n",
    "                if val < 0:\n",
    "                    sign = -1\n",
    "                    \n",
    "                if col == i or splitting[col] == 1:\n",
    "                    continue\n",
    "                coarse_sum = 0\n",
    "                for k in range(A.indptr[col], A.indptr[col+1]):\n",
    "                    colk = A.indices[k]\n",
    "                    if pos[colk] >= 0 or colk == i:\n",
    "                        val = A.data[k]\n",
    "                        if val * sign < 0:\n",
    "                            coarse_sum += val\n",
    "                if coarse_sum == 0:\n",
    "                    weak_sum += S.data[j]\n",
    "                else:\n",
    "                    coarse_sum = row_strong[col] / coarse_sum\n",
    "                    for k in range(A.indptr[col], A.indptr[col+1]):\n",
    "                        colk = A.indices[k]\n",
    "                        idx = pos[colk]\n",
    "                        if idx >= 0:\n",
    "                            val = A.data[k]\n",
    "                            if val * sign < 0:\n",
    "                                Px[idx] += coarse_sum * val \n",
    "                        if colk == i:\n",
    "                            weak_sum += coarse_sum * val\n",
    "            \n",
    "            for j in range(S.indptr[i], S.indptr[i+1]):\n",
    "                col = S.indices[j]\n",
    "                row_strong[col] = 0\n",
    "                                \n",
    "            for j in range(Pp[i], Pp[i+1]):\n",
    "                Px[j] /= -weak_sum\n",
    "                col = Pj[j]\n",
    "                pos[col] = -1\n",
    "\n",
    "                            \n",
    "    cols = np.zeros(nnz, dtype = 'int')\n",
    "    for i in range(nnz):\n",
    "        cols[i] = Pj[i]\n",
    "    cols.sort()\n",
    "    \n",
    "    cols_to_new = np.empty(A.shape[0], dtype='int')\n",
    "    cols_to_new.fill(-1)\n",
    "    \n",
    "    n_cols = 0\n",
    "    col = cols[0]\n",
    "    cols_to_new[col] = n_cols\n",
    "    n_cols += 1\n",
    "    for i in range(1, nnz):\n",
    "        col = cols[i]\n",
    "        if col != cols[n_cols-1]:\n",
    "            cols_to_new[col] = n_cols\n",
    "            cols[n_cols] = col\n",
    "            n_cols += 1\n",
    "            \n",
    "    for i in range(nnz):\n",
    "        Pj[i] = cols_to_new[Pj[i]]\n",
    "                    \n",
    "    return csr_matrix((Px, Pj, Pp), shape = (A.shape[0], n_cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MIS-2 Algorithm\n",
    "def mis2(A):\n",
    "    n = A.shape[0]\n",
    "    \n",
    "    # 1\n",
    "    np.random.seed(100) # For comparisons\n",
    "    r = np.random.rand(n)\n",
    "    \n",
    "    #2\n",
    "    A_dir = A.copy().tocsr()\n",
    "    for i in range(n):\n",
    "        for j in range(A_dir.indptr[i], A_dir.indptr[i+1]):\n",
    "            col = A_dir.indices[j]\n",
    "            if (r[i] <= r[col]):\n",
    "                A_dir.data[j] = 0.0\n",
    "    A_dir.eliminate_zeros()\n",
    "\n",
    "    # 3\n",
    "    S = []\n",
    "    \n",
    "    # 4\n",
    "    V = np.arange(n)\n",
    "    states = np.full((n,1), -1)\n",
    "    remaining = n\n",
    "        \n",
    "    iterate = 0\n",
    "    \n",
    "    while remaining > 0:\n",
    "        # 6\n",
    "        print(\"Iteration\", iterate, \"# active nodes\", len(V), \"/\", n, \", MIS2 set size is\", len(S))\n",
    "        iterate+=1\n",
    "        L = []\n",
    "        for v in V:\n",
    "            found = False\n",
    "            for jj in range(A_dir.indptr[v], A_dir.indptr[v+1]):\n",
    "                w = A_dir.indices[jj]\n",
    "                if states[w] < 0:\n",
    "                    found = True\n",
    "                    break\n",
    "            if not found:\n",
    "                L.append(v)\n",
    "                states[v] = -2\n",
    "                        \n",
    "        # 8\n",
    "        for v in L:\n",
    "            found = False\n",
    "            for jj in range(A.indptr[v], A.indptr[v+1]):\n",
    "                w = A.indices[jj]\n",
    "                for kk in range(A.indptr[w], A.indptr[w+1]):\n",
    "                    u = A.indices[kk]\n",
    "                    if states[u] < -1 and r[u] > r[v]:\n",
    "                        found = True\n",
    "                        break\n",
    "            if not found:\n",
    "                S.append(v)\n",
    "                states[v] = -3                    \n",
    "        \n",
    "        # 11\n",
    "        for v in V:\n",
    "            if states[v] == -3:\n",
    "                continue\n",
    "            found = False\n",
    "            for jj in range(A.indptr[v], A.indptr[v+1]):\n",
    "                u = A.indices[jj]\n",
    "                if states[u] == -3:\n",
    "                    found = True\n",
    "                    break\n",
    "            if not found:\n",
    "                for jj in range(A.indptr[v], A.indptr[v+1]):\n",
    "                    w = A.indices[jj]\n",
    "                    for kk in range(A.indptr[w], A.indptr[w+1]):\n",
    "                        u = A.indices[kk]\n",
    "                        if states[u] == -3:\n",
    "                            found = True\n",
    "                            break\n",
    "                    if found:\n",
    "                        break\n",
    "            if found:\n",
    "                states[v] = -4\n",
    "            \n",
    "        ctr = 0\n",
    "        for i in range(remaining):\n",
    "            v = V[i]\n",
    "            if (states[v] == -4):\n",
    "                states[v] = 0\n",
    "            elif (states[v] == -3):\n",
    "                states[v] = 1\n",
    "            else:\n",
    "            #if (states[V[i]] < 0):\n",
    "                V[ctr] = V[i]\n",
    "                ctr += 1\n",
    "        V.resize(ctr)\n",
    "        remaining = ctr\n",
    "        \n",
    "    print(\"Final MIS2 set size is\", len(S))\n",
    "\n",
    "    return S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "###    Gallery Tests\n",
    "###############################\n",
    "# Rotated Anisotropic Diffusion\n",
    "n = 25\n",
    "grid = (n,n)\n",
    "sten = diffusion_stencil_2d(epsilon=0.001, theta=np.pi/8)\n",
    "Aniso = stencil_grid(sten, grid, dtype='float').tocsr()\n",
    "PetscPrint(\"aniso.pm\", Aniso)\n",
    "\n",
    "# Laplacian 27pt\n",
    "n = 10\n",
    "grid = (n,n,n)\n",
    "sten = [[[-1,-1,-1],[-1,-1,-1],[-1,-1,-1]],\n",
    "       [[-1,-1,-1],[-1,26,-1],[-1,-1,-1]],\n",
    "       [[-1,-1,-1],[-1,-1,-1],[-1,-1,-1]]]\n",
    "Laplacian = stencil_grid(sten, grid, dtype='float').tocsr()\n",
    "PetscPrint(\"laplacian.pm\", Laplacian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "###    Core Tests\n",
    "###############################\n",
    "# Transposes\n",
    "AT = csr_matrix(Aniso)\n",
    "AT.T\n",
    "PetscPrint(\"aniso_T.pm\", AT)\n",
    "\n",
    "AT = csr_matrix(Laplacian)\n",
    "AT.T\n",
    "PetscPrint(\"laplacian_T.pm\", AT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "###    Strength Tests\n",
    "###############################\n",
    "Aniso_S = ClassicStrength(Aniso, 0.25)\n",
    "PetscPrint(\"aniso_S.pm\", Aniso_S)\n",
    "\n",
    "Laplacian_S = ClassicStrength(Laplacian, 0.25)\n",
    "PetscPrint(\"laplacian_S.pm\", Laplacian_S)\n",
    "\n",
    "Aniso_SS = SymmetricStrength(Aniso, 0.25)\n",
    "PetscPrint(\"aniso_SS.pm\", Aniso_SS)\n",
    "\n",
    "Laplacian_SS = SymmetricStrength(Laplacian, 0.25)\n",
    "PetscPrint(\"laplacian_SS.pm\", Laplacian_SS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 125) 8792\n"
     ]
    }
   ],
   "source": [
    "################################\n",
    "###    Ruge Stuben AMG Tests\n",
    "###############################  \n",
    "# RS Splitting\n",
    "AnisoSplit = RS(Aniso_S)\n",
    "np.savetxt(\"aniso_split.txt\", AnisoSplit, fmt=\"%d\")\n",
    "LaplacianSplit = RS(Laplacian_S)\n",
    "np.savetxt(\"laplacian_split.txt\", LaplacianSplit, fmt=\"%d\")\n",
    "\n",
    "# Direct Interpolation\n",
    "from pyamg.classical.interpolate import direct_interpolation\n",
    "Aniso_P = direct_interpolation(Aniso, Aniso_S, AnisoSplit)\n",
    "PetscPrint(\"aniso_P_direct.pm\", Aniso_P)\n",
    "Laplacian_P = direct_interpolation(Laplacian, Laplacian_S, LaplacianSplit)\n",
    "PetscPrint(\"laplacian_P_direct.pm\", Laplacian_P)\n",
    "\n",
    "# Modified Classical Interpolation\n",
    "Aniso_P = ModClassicalInterp(Aniso, Aniso_S, AnisoSplit)\n",
    "PetscPrint(\"aniso_P_mod_class.pm\", Aniso_P)\n",
    "Laplacian_P = ModClassicalInterp(Laplacian, Laplacian_S, LaplacianSplit)\n",
    "PetscPrint(\"laplacian_P_mod_class.pm\", Laplacian_P)\n",
    "\n",
    "# Extended+i Interpolation\n",
    "Aniso_P = ExtendInterpolation(Aniso, Aniso_S, AnisoSplit)\n",
    "PetscPrint(\"aniso_P_extend.pm\", Aniso_P)\n",
    "Laplacian_P = ExtendInterpolation(Laplacian, Laplacian_S, LaplacianSplit)\n",
    "PetscPrint(\"laplacian_P_extend.pm\", Laplacian_P)\n",
    "print(Laplacian_P.shape, Laplacian_P.nnz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 # active nodes 625 / 625 , MIS2 set size is 0\n",
      "Iteration 1 # active nodes 152 / 625 , MIS2 set size is 31\n",
      "Iteration 2 # active nodes 59 / 625 , MIS2 set size is 44\n",
      "Iteration 3 # active nodes 17 / 625 , MIS2 set size is 50\n",
      "Iteration 4 # active nodes 8 / 625 , MIS2 set size is 51\n",
      "Iteration 5 # active nodes 2 / 625 , MIS2 set size is 52\n",
      "Final MIS2 set size is 53\n",
      "Iteration 0 # active nodes 1000 / 1000 , MIS2 set size is 0\n",
      "Iteration 1 # active nodes 109 / 1000 , MIS2 set size is 19\n",
      "Iteration 2 # active nodes 31 / 1000 , MIS2 set size is 27\n",
      "Iteration 3 # active nodes 1 / 1000 , MIS2 set size is 33\n",
      "Final MIS2 set size is 34\n"
     ]
    }
   ],
   "source": [
    "################################\n",
    "###    Smoothed Aggregation Tests\n",
    "###############################        \n",
    "# Test MIS\n",
    "Aniso = Aniso.tocsr()\n",
    "S = mis2(Aniso)\n",
    "S = np.array(S)\n",
    "split_Aniso = np.zeros(Aniso.shape[0])\n",
    "for s in S:\n",
    "    split_Aniso[s] = 1\n",
    "np.savetxt(\"aniso_mis.txt\", split_Aniso, fmt=\"%d\")\n",
    "\n",
    "Laplacian = Laplacian.tocsr()\n",
    "S = mis2(Laplacian)\n",
    "S = np.array(S)\n",
    "split_Laplacian =  np.zeros(Laplacian.shape[0])\n",
    "for s in S:\n",
    "    split_Laplacian[s] = 1\n",
    "np.savetxt(\"laplacian_mis.txt\", split_Laplacian, fmt=\"%d\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
